home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / fit2d.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  3KB  |  130 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  *    fit2d - 
  19.  *        fit a plane to a 2D scalar function.
  20.  *
  21.  *                Paul Haeberli with Tom Davis
  22.  */
  23. #include "stdio.h"
  24.  
  25. static double m00, m01, m02;
  26. static double m10, m11, m12;
  27. static double m20, m21, m22;
  28. static double v0, v1, v2;
  29. static double A, B, C;
  30. static int havemat;
  31.  
  32. initfit()
  33. {
  34.     m00 = m01 = m02 = 0.0;
  35.     m10 = m11 = m12 = 0.0;
  36.     m20 = m21 = m22 = 0.0;
  37.     v0 = v1 = v2 = 0.0;
  38.     havemat = 0;
  39. }
  40.  
  41. fitpoint(x,y,r)
  42. float x, y, r;
  43. {
  44.     m00 += x*x;
  45.     m01 += x*y;
  46.     m02 += x;
  47.     m10 += x*y;
  48.     m11 += y*y;
  49.     m12 += y;
  50.     m20 += x;
  51.     m21 += y;
  52.     m22 += 1.0;
  53.  
  54.     v0 += r*x;
  55.     v1 += r*y;
  56.     v2 += r;
  57.     havemat = 0;
  58. }
  59.  
  60. float getfitval(x,y)
  61. float x, y;
  62. {
  63.     float inpos, outpos; 
  64.     float mat[4][4], imat[4][4];
  65.     float a, b, c;
  66.  
  67.     if(!havemat) {
  68.     bzero(mat,16*sizeof(float));
  69.     mat[0][0] = m00;
  70.     mat[0][1] = m01;
  71.     mat[0][2] = m02;
  72.     mat[1][0] = m10;
  73.     mat[1][1] = m11;
  74.     mat[1][2] = m12;
  75.     mat[2][0] = m20;
  76.     mat[2][1] = m21;
  77.     mat[2][2] = m22;
  78.       mat[3][3] = 1.0;
  79.     if(m22<0.5) {
  80.         A = 0.0;
  81.         B = 0.0;
  82.         C = 0.0;
  83.     } else if(invertmat(mat,imat)) {
  84.         xformpnt(imat,v0,v1,v2,&a,&b,&c);
  85.         A = a;
  86.         B = b;
  87.         C = c;
  88.     } else {
  89.         A = 0.0;
  90.         B = 0.0;
  91.         C = v2/m22;
  92.     }
  93.     havemat = 1;
  94.     }
  95.     return A*x+B*y+C;
  96. }
  97.  
  98. #ifdef TESTIT
  99. printval(x,y)
  100. float x, y;
  101. {
  102.     printf("val iat %f %f is %f\n",x,y,getfitval(x,y));
  103. }
  104.  
  105. main()
  106. {
  107.     initfit();
  108.     fitpoint(0.0,0.0,0.0);
  109.     fitpoint(1.0,0.0,1.0);
  110.     fitpoint(1.0,0.0,1.0);
  111.     fitpoint(1.0,0.0,1.0);
  112.     fitpoint(1.0,0.0,1.0);
  113.     fitpoint(1.0,0.0,1.0);
  114.     fitpoint(1.0,0.0,1.0);
  115.     fitpoint(1.0,0.0,1.0);
  116.     fitpoint(1.0,0.0,1.0);
  117.     fitpoint(1.0,0.0,1.0);
  118.     fitpoint(1.0,0.0,1.0);
  119.     fitpoint(1.0,0.0,1.0);
  120.     fitpoint(0.0,1.0,1.0);
  121.     fitpoint(1.0,1.0,2.0);
  122.     printval(0.0,0.0);
  123.     printval(1.0,0.0);
  124.     printval(0.0,1.0);
  125.     printval(1.0,1.0);
  126.     printval(0.1,0.0);
  127.     printval(0.5,0.5);
  128. }
  129. #endif
  130.